//----------------------------------------------------------------------------
// M A R I T I M E  R E S E A R C H  I N S T I T U T E  N E T H E R L A N D S
//----------------------------------------------------------------------------
// Copyright (C) 2008 - MARIN - All rights reserved - http://www.marin.nl
//----------------------------------------------------------------------------
// Program    : mermaid
// Module     : vbm
// File       : WeakRefl.cpp
//----------------------------------------------------------------------------

#include <algorithm> // min,max
#include "FloatUtils.h"
#include "WeakRefl.h"
#include "assertplus.h"

#ifdef _WINDOWS
#pragma warning(disable: 4305) // truncation from 'double' to 'float'
#endif

/// constructor
WeakRefl::WeakRefl(int somfeld, int Nx1, int Nx2)
: m_somfeld(somfeld)
, m_mu_num(5.0 / 6.0)
, m_mu_den(1.0 / 3.0)
{
}

/// destructor
WeakRefl::~WeakRefl()
{
}

/// computes weak reflections
void WeakRefl::compute(Array2D<REAL> &zetaN,
                       Array2D<REAL> &QzetaN,
                       const Array1D<REAL> &zetaC,
                       const Array1D<REAL> &zetaS,
                       const Array1D<REAL> &cN,
                       int Nx,
                       REAL dx2,
                       REAL dx1,
                       REAL dt)
{
    // shift the old values and
    // update Q-values with Sommerfeld boundary condition
    for (int i = 1; i <= Nx; ++i)
    {
#ifdef indices
        zetaN [i][2] = zetaN [i][1];
        zetaN [i][1] = zetaN [i][0];
        QzetaN[i][2] = QzetaN[i][1];
        QzetaN[i][1] = QzetaN[i][0];
        QzetaN[i][0] = QzetaN[i][1] - cN[i] * dt / (2 * dx2) * (3 * zetaN[i][1] - 4 * zetaC[i] + zetaS[i]);
#else
        REAL* pzeta  = zetaN [i];
        REAL* pqzeta = QzetaN[i];
        *(pzeta + 2) = *(pzeta + 1);
        *(pzeta + 1) = *pzeta;
        *(pqzeta + 2) = *(pqzeta + 1);
        *(pqzeta + 1) = *pqzeta;
        *pqzeta = *(pqzeta + 1) - cN[i] * dt / (2 * dx2) * (3 * *(pzeta + 1) - 4 * zetaC[i] + zetaS[i]);
#endif
    }
    QzetaN[0     ][0] = (4 * QzetaN[1 ][0] - QzetaN[2     ][0]) / 3;
    QzetaN[Nx + 1][0] = (4 * QzetaN[Nx][0] - QzetaN[Nx - 1][0]) / 3;

    // update the boundary values
    switch (m_somfeld)
    {

        case 1: // Classical Sommerfeld condition
            for (int i = 1; i <= Nx; ++i)
            {
                zetaN[i][0] = QzetaN[i][0];
            }
            break;

        case 2: // Sommerfeld with improved directional characteristics
        {
            Array1D<REAL> tmp(Nx + 2);
            Array1D<REAL> curv(Nx + 1);
            // NOTE ASR: this leads to very strange results, in the original Fortran code too!
            for (int i = 0; i <= Nx + 1; ++i)
            {
                tmp[i] = m_mu_num * zetaN[i][1] - m_mu_den * QzetaN[i][0];
            }
            for (int i = 1; i <= Nx; ++i)
            {
                curv[i] = tmp[i - 1] - 2 * tmp[i] + tmp[i + 1];
            }
            curv[1] = 2 * tmp[1] - 5 * tmp[2] + 4 * tmp[3] - tmp[4];
            curv[Nx] = 2 * tmp[Nx] - 5 * tmp[Nx - 1] + 4 * tmp[Nx - 2] - tmp[Nx - 3];
            for (int i = 1; i <= Nx; ++i)
            {
                zetaN[i][0] = 2 * zetaN[i][1] - zetaN[i][2] +
                    (QzetaN[i][0] - 2 * QzetaN[i][1] + QzetaN[i][2]) +
                    FloatUtils::square(cN[i] * dt / dx1) * curv[i];
            }
            //  classical Sommerfeld in the corners
            zetaN[0][0] = (4 * zetaN[1][0] - zetaN[2 ][0]) / 3;
            // Neumann b.c.
            zetaN[Nx + 1][0] = (4 * zetaN[Nx][0] - zetaN[Nx - 1][0]) / 3;
            break;
        }
        default:
            ASSERTFALSE();
    }
}
